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We investigate the sandpile model on the two-dimensional Sierpinski gasket fractal. We find that 
the model displays novel critical behavior, and we analyze the distribution functions of avalanche 
sizes, lifetimes and topplings and calculate the associated critical exponents r = 1.51 ± 0.04, ct = 
1.63±0.04 and /i = 1.36±0.04. The avalanche size distribution shows power law behavior modulated 
by logarithmic oscillations which can be related to the discrete scale invariance of the underlying 
lattice. Such a distribution can be formally described by introducing a complex scaling exponent 
T* = r + j(5, where the real part r corresponds to the power law and the imaginary part 5 is related 
to the period of the logarithmic oscillations. 
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I. INTRODUCTION 

The concept of self-organized criticality (SOC) has 
been introduced by Bak et al. to describe the ten- 
dency of a large class of dynamical systems to sponta- 
neously evolve into a critical state without fine tuning 
of any external parameter. Sandpile models have 
been introduced as an example of this kind of phenom- 
ena and have been widely studied numerically and ana- 
lytically . Two principal analytical approaches have 
been followed: the first involves the group theory formal- 
ism introduced by Dhar |^ and the second is a real space 
renormalization scheme recently developed by Pietronero 
et al. . Other theoretical approaches involve nonlinear 
continuous differential equations pO| , pT[ |. 

Sandpile models are inspired by the dynamics of sand 
flowing along the slope of a pile. By adding sand grains to 
the pile the system eventually reaches a stationary state 
characterized by avalanches of all length scales. The term 
criticality refers here to the absence of any characteris- 
tic length scale in this state. Sandpile models have been 
studied mostly on Euclidean lattices. It has been shown 
that different kinds of Euclidean lattices do not affect the 
critical exponents This fact is similar to the univer- 
sality observed in ordinary critical phenomena. More- 
over, in the case of the Bethe lattice one recovers the 
mean field results [ p^lJT^ . However, sandpile models, to 
our knowledge, have not been studied on a fractal sub- 
strate, in particular on a simple deterministic fractal such 
as that epitomized by the Sierpinski gasket (SG). 

It has been shown, via specific calculations and 
through general rigorous analysis , that for the stan- 
dard Ising model (and for some more general models) on 
finitely ramified fractals no spontaneous magnetization 
can exist at any finite temperature. It might have hap- 
pened that, by some assumed analogy, no self-organized 
critical behavior has been expected so far to occur on the 
deterministic fractals. However, we shall demonstrate 
that the SOC phenomenon exists on the SG fractal and 
displays novel features. Specifically, we study numeri- 



cally the critical height sandpile model on the SG lattice 
with the generator scaling base 6 = 2 which corresponds 
to the fractal dimension D = log 3/ log 2 « 1.58. We 
calculate the distribution of avalanche sizes, their life- 
times, and topplings. The avalanche size distribution 
shows a power law behavior modulated by logarithmic 
oscillations. This kind of oscillation has been already ob- 
served in the scaling functions of different systems , 
and here it can be related to the discrete scale invariant 
nature of the underlying fractal lattice. It is interesting to 
note that complex scaling exponents have been recently 
detected in earthquakes statistics . 

The measured scaling exponents vary with the system 
size L and the values, extrapolated to L — > oo, appear 
to differ from those computed on the Euclidean lattices. 
Computing expectation values, we are able to verify the 
relationships between different critical exponents. 

In addition, we investigate time correlations of the 
number of drops and topplings during the avalanche. 
Calculating the power spectra, we find that as in the case 
of the two-dimensional Euclidean lattice there are 

no long-range temporal correlations. 



II. THE MODEL 

Our cellular automata model is defined on the SG lat- 
tice as shown in Fig. 1. The number n is related to the 
number of sites L = 2" + 1 along one direction of the 
lattice and is used hereafter as a measure of the system 
size. Within the sandpile model, all the sites of the frac- 
tal lattice are exposed to the same local dynamical rules. 
The exceptions are the three apex sites where the sand 
grains flow out of the system. The dynamics begins when 
we associate an integer height variable Zi with every lat- 
tice site i. At each later step one lattice site is chosen 
at random and its height is increased by one. Whenever 
the height on a site i reaches the critical value = 4, 
the site becomes unstable (active) and relaxes according 
to the following rules: 
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Zi- 4, (1) 
Zj — > + 1, (2) 

where j are the nearest neighbors of i. These rules con- 
serve the total number of grains, except on the three 
apex sites (independent of the system size) where two 
sand grains are lost. Successive relaxation events gener- 
ate the sand flow that eventually brings the sand out of 
the system. Due to the local conservation, imposed by 
the dynamical rules, the system finally evolves into the 
stationary state characterized by the balance between the 
input and the output flow. 

The critical exponents are extracted from avalanche 
distributions in the stationary state. We define the size 
S as the number of distinct sites visited by an avalanche, 
the toppling size m as the number of relaxation events 
and the lifetime T as the number of updating steps dur- 
ing an avalanche. All these quantities are expected to be 
distributed as power laws 

P{S) ^ S-\ (3) 
P(m) - m^'', (4) 
P{T) - (5) 

where r, fi, and a, are critical exponents of the respective 
distribution functions P{S), P{m) and P{T). We can 
relate these exponents by considering conditional expec- 
tation values for an average avalanche size < St > and 
an average toppling size < ttit > at a given avalanche 
lifetime T: 

< St >^ , (6) 

< TUT >~ , (7) 

and similar other relations. By taking into account the 
definitions of critical exponents, given by Eqs.(3)-(7), 
scaling relations between exponents can be derived ||6|]: 

r = l + (a-l)//3, (8) 

M=l + (a-l)/7. (9) 

Finally we study temporal correlations by consider- 
ing the number of active sites and the number of grains 
falling out of the system at each time step. The power 
spectrum of this signal falls off as 

S{f) ~ /-^ (10) 

In the Euclidean case, = 2 showing the absence of 
long-range temporal correlations. 

III. SIMULATION RESULTS 

We perform numerical simulations for different lattice 
sizes ranging from n = 3 to n = 7. The total number 



of sites Sn+i of the system size n -f 1 is related to the 
number of sites Sn of the system size n via the equation 
Sn+i = SS*.,! — 3 with So — 3, which corresponds to a total 
number of lattice sites going from S3 = 42 to Sj = 3282. 

A simple way to characterize the properties of the sta- 
tionary state is to compute the fraction of sites having 
height Zi = 1, 2, 3 and 4. We report these results in Ta- 
ble 1 for different system sizes, together with an average 
height < z >. The obtained values are very close to those 
found on the Euclidean lattice |^,|^ . 

In Fig. 2 we show the avalanche size distribution for 
different system sizes. One can see that there are quite 
a few avalanches (represented by the last peak of each 
distribution curve) which span the entire lattice. This 
occurs because the balance between incoming and out- 
going particles forces the avalanches to reach the three 
apex sites. Due to the self-similarity of the underlying 
lattice, the same phenomenon occurs on all fractal sub- 
structures, which is manifested by a series of peaks on 
each distribution curve. 

The phenomenon described above is reflected in a pe- 
culiar behavior of the active sites during the evolution 
of an avalanche: the active sites are localized (trapped) 
within fractal substructures for many time steps. Such 
a trapping does not occur in the Euclidean lattice where 
the active sites are essentially on the avalanche front. 
This is apparent from Fig. 3 where the active sites for a 
typical avalanche in the Euclidean square lattice are com- 
pared with an avalanche on the fractal lattice. Similar 
differences which spring from differences in the topology 
of the lattices were noted before ||2^ in a study of linear 
polymers on the diamond hierarchical lattice. 

The power-law behavior in a double logarithmic plot 
is modulated by oscillations with a period p that can be 
related to the scaling properties of the SG lattice. A 
self-similar lattice is left invariant only by a discrete set 
of scale transformations, namely by those with a scaling 
parameter of the form A = 6". Under this condition, it 
has been shown |2^ that the most general scale invari- 
ant function of the real space coordinates is a power law 
multiplied by a logarithmically periodic function. These 
oscillations can be formally described by introducing a 
complex scaling exponent t* = t + iS where the real 
part T corresponds to the power law exponent, while the 
imaginary part 6 is related to the period of oscillations. 
In our case S = 27r/p = 27r/log3 = 5.72. To extract r, 
we fit the distribution with a power law modulated by a 
periodic function. 

The last peak in the distribution of avalanche sizes is a 
consequence of the fact that at any system size there are 
only three boundary points where the sand can flow out 
of the system, in contrast to the Euclidean lattice where 
the number of boundary points increase in proportion to 
the system size. However, one can study the effect of 
boundaries by calculating the avalanche size distribution 
of the same sandpile model on the Euclidean lattice with 
only four boundary points, e.g. the four corner sites of 
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the square lattice (for the rest of the edge points peri- 
odic boundary conditions apply). The results are shown 
in Fig. 4 for system sizes L — 4, 8, 16 and 32. The distri- 
butions are power laws with peaks at the total numbers 
of sites on the lattice. 

In the SG case, we report in Figs. 5 and 6 the distri- 
butions of avalanche lifetimes and topplings in a double 
logarithmic plot. Both distributions display pure power- 
law behavior without any modulations. The power-law 
regimes grow with the system size. 

As in the case of the Euclidean lattice Q , the scaling 
exponents depend on the system size. We can, however, 
extract the asymptotic results by plotting the logarithm 
of the exponents versus 1 / log L, where L is the linear size 
of the lattice. This relationship is presented in Fig. 7, 
where we depict also the extrapolated critical exponents 
in the limit n — > cx). We found the following extrapo- 
lated values Too = 1-51 ± 0.04, Uoo = 1-64 ± 0.04 and 
/ioo = 1.36 ±0.04 for the avalanche size, lifetime and top- 
pling distributions, respectively. The distribution func- 
tions presented in these figures have been calculated by 
averaging over 2^^ avalanches. 

To check the consistency of our results we compute 
the scaling exponents of the conditional expectation val- 
ues, i.e. exponents related to the average avalanche size 
< St > and number of topplings < jtlt > in dependence 
on the lifetime T. Fig. 8 shows results of our compu- 
tation in a double logarithmic plot. The slopes in the 
figure correspond to the exponents (3 and 7 as defined 
in Eqs.(|l) and @) and are found to be ^ = 1.13 ± 0.05 
and 7 — 1.73 ± 0.05. These two values can be com- 
pared with the ones evaluated from the scaling relations 
given by Eqs.(8) and (9), and using the estimated crit- 
ical exponents Too, aoo and ^00 • Thus, we find the val- 
ues P = 1.24 ± 0.12 and 7 = 1.75 ± 0.18, which are in 
agreement, within the numerical error, with the directly 
obtained values from Fig. 8. 

Finally, within the scope of the sandpile model ||l[ , we 
calculate the temporal correlations of two quantities: the 
number of particles which fall out of a system in a unit 
time and the number of topplings. The unit time in this 
case corresponds to one updating step of the lattice vari- 
ables. We calculate the power spectra of the above two 
quantities. In both cases we find a 1//^ type of spec- 
trum. Thereby, the type of temporal correlations turns 
out to be the same as in the original model on a two- 
dimensional square lattice [|l8|,^. Our results are pre- 
sented in Fig. 9. The flattening of the 1//'^ spectrum at 
small frequencies is due to finite size effects. In contrast 
to other scaling exponents presented in this paper, the 
exponents of power spectra do not vary with the system 
size, that is, the 1//^ type of spectrum corresponds to an 
exponential decay of temporal correlations independently 
of the system size. 



IV. CONCLUSION 

We computed numerically the scaling exponents for 
the avalanche distributions in the critical height sandpile 
model on the Sierpinski gasket lattice with the generator 
base b — 2. The lattice coordination number is the same 
as in the two-dimensional square lattice and therefore 
the dynamical rules of the sandpile model are exactly 
the same. The boundaries, however, are different, since 
on the SG lattice the sand can flow out only through 
three sites at every scale. This fact changes substan- 
tially the avalanche dynamics. The active sites become 
trapped (localized) and topple more than once during a 
single avalanche. 

In relation to the standard critical phenomena it is in- 
teresting to note how the dimensionality and the topology 
of the lattice affect the critical behavior of the model. In 
one dimension, the critical height sandpile model is triv- 
ial in that avalanches are not power-law distributed [Q. 
A similar behavior occurs for example in the Ising model 
where no phase transition is observed in dimension less 
than two. We have shown, however, that on a finitely 
ramified fractal, the sandpile has nontrivial critical be- 
havior, in contrast to the Ising model which has no phase 
transition on such a fractal lattice Qjl5| . 

Finally, we note that self-similar lattices have been 
proven very helpful in constructing exact real space 
renormalization group transformations p^j2^ for stan- 
dard critical phenomena. Having demonstrated that self- 
organized criticality can exist on a fractal lattice, it would 
be beneficial to find such a transformation for this model, 
trying to link the rigorous approach of Dhar et al. j^] 
with the real-space renormalization scheme presented in 
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FIG. 4. The distributions of avalanche sizes for a sandpile 
model on a square lattice with four exiting points are depicted 
for several system sizes L = 4, 8, 16 and 32. Data for L = 16 
and L — 32 are averaged at sizes S > 50. 



FIG. 5. The distribution of avalanche lifetimes of the sand- 
pile model on a SG lattice as calculated for different system 
sizes. Data are logarithmically binned at lifetimes T > 50. 



FIG. 6. The distribution of the number of topplings per 
avalanche of the sandpile model on a SG lattice as calculated 
for different system sizes. Data are logarithmically binned at 
toppling sizes m > 50. 



FIG. 7. The logarithm of the critical exponents found as a 
result of a simulation for different lattice sizes plotted against 
1/n. Estimations of the three exponents in the limit of in- 
finitely large lattice size (1/n — > 0) are also shown. The ex- 
trapolated exponents in the limit n oo are Too = 1.51±0.04, 
Qfoc = 1.63 ± 0.04 and /loo = 1.36 ± 0.04. 



FIG. 8. The average avalanche size < St > and the aver- 
age number of topplings < ttit > as functions of the lifetime 
T, presented in a double logarithmic plot. The corresponding 
slopes are /3 = 1.13 ± 0.05 and 7 = 1.75 ± 0.05. Data are 
logarithmically binned at lifetimes T > 50. 



FIG. 9. Power spectra of the number of particles which 
drop out of the system (stars) and of the number of topplings 
(diamonds) in a given time step. The data for three system 
sizes are depicted, from n = 5 to n = 7. The results show 
l/f^ type of spectra, which indicates an exponential decay 
of time correlations. It can be observed that the frequency 
at which the 1//^ type of the spectrum crosses over to the 
white noise type (i.e. the flat part of the spectrum), due to 
the finite size of the system, decreases with the lattice size. 



FIG. 1. An example of the SG lattice with the generator 
6 = 2, at the stage of construction n = 2, and with linear size 
L = 5. Each site has four neighbors except for the three apex 
sites with only two neighboring sites. The arrows indicate the 
direction of sand flow from a chosen site. 

FIG. 2. The distribution of avalanche sizes of the sand- 
pile model on a SG lattice. Different curves correspond to 
different system sizes. The arrows indicate the peaks in the 
distribution. 
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TABLE I. The fraction of sites having height equal to 
z = 1, 2, 3 and 4 and the average height (z) as found for 
different system sizes L = 2" -f 1 and n = 3, 4, 5, 6 and 7. 



FIG. 3. A snapshot of an avalanche on the Sierpinski gas- 
ket lattice (a) compared with an avalanche on the Euclidean 
square lattice (b). Active sites are depicted in black and sites 
that have toppled at least once are colored in gray. 
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